
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here
    clear;
    clc;
    setEnv();
    workingDir = 'C:\Users\rniehaus\Documents\Elegans\';
    studyInstanceName = 'Study-RFUM294_loops';
    [dataList, suppData] = getStudyInstance(workingDir, studyInstanceName);

    for iFrame = 1:length(dataList)
        bwImage = suppData(iFrame).BwImage;
        
          %Plot bw image
             figure(1), imshow(bwImage) 
        
        %% GEt outer boundary
        [B,L,N,A] = bwboundaries(bwImage);

        outerBoundary = B{1};

        %Convert from subscripts to coordinates
        vidRows = size(bwImage,1);
        yOuter = vidRows - outerBoundary(:,1);
        xOuter = outerBoundary(:,2);
        outer = [xOuter,yOuter];
        
        %% Get inner boundary
            
        if length(B) > 1;
            innerBoundaryCell = B(2);
            innerBoundaryMat = innerBoundaryCell{1};
            yInner = vidRows - innerBoundaryMat(:,1);
            xInner = innerBoundaryMat(:,2);
            inner = [xInner, yInner];
             totalBoundary = [[xOuter;xInner],[yOuter; yInner] ];
        else
            xInner = [];
            yInner = [];
            inner = [];
            totalBoundary = [xOuter,yOuter ];
        end
        
       
            
        %% Delaunay

        tri = delaunay(totalBoundary(:,1), totalBoundary(:,2));

        %Plot the entire boundary and mesh
        figure(3)
        triplot(tri, totalBoundary(:,1), totalBoundary(:,2)); 
            
        
        %% Get side 1 and side 2

        %Get Sides A (from head to tail) and Side B(from head to tail)of
        %the outer boundary
        
        
        if dataList(iFrame).Curv_EndARow == 1 && dataList(iFrame).Curv_EndACol == 1
            %Delta configuration
            endA = [dataList(iFrame).Curv_EndBCol, vidRows - dataList(iFrame).Curv_EndBRow ];
            dist = round(size(outer,1)/2);
            endB = [outer(dist,1), outer(dist,2)];
        elseif dataList(iFrame).Curv_EndBRow == 1 && dataList(iFrame).Curv_EndBCol == 1
            %Delta configuration
            endA = [dataList(iFrame).Curv_EndACol, vidRows - dataList(iFrame).Curv_EndARow ];
            dist = round(size(outer,1)/2);
            endB = [outer(dist,1), outer(dist,2)];
        else    
          
            endA = [dataList(iFrame).Curv_EndACol, vidRows - dataList(iFrame).Curv_EndARow ];
            endB = [dataList(iFrame).Curv_EndBCol, vidRows - dataList(iFrame).Curv_EndBRow ];
        end 
        
        %Find locations in boundary
        [~, locA] = ismember(endA, outer, 'rows');
        [~, locB] = ismember(endB, outer, 'rows');
        
        side1 = [];
        side2 = [];

        if locA < locB
            side1 = outer(locA:locB,:);
            side2 = [outer(locB:size(outer,1),:); outer(1:locA,:)];
        else
            side1 = outer(locB:locA,:);
            side2 = [outer(locA:size(outer,1),:); outer(1:locB,:)];
        end
       
            

        %Plot inner and outer boundary with sides shown
        figure(4)
        plot(side1(:,1), side1(:,2), 'g');
        hold on
        plot(side2(:,1), side2(:,2), 'm');
        hold on
        if ~isempty(xInner)
            plot(xInner, yInner, 'b'); 
        end
        hold off
            
       %List of boundary points referenced    
       refBoundaryPts = unique(reshape(tri,size(tri,1)*3, 1)  ) ;
       outerOnly = refBoundaryPts(refBoundaryPts <= size(outer,1) );
       isReferenced = zeros(size(totalBoundary,1),1);
       isReferenced( refBoundaryPts) = 1;
       
       
%       isMember
%        for i = 1:size(refPtsOnOuter,1)
%            col1 = find(tri(:,1)== refPtsOnOuter(i));
%            col2 = find(tri(:,2)== refPtsOnOuter(i));
%            col3 = find(tri(:,3)== refPtsOnOuter(i));
%        end
%        total = col1 + col2 + col3;
%        hasOuterEdges  = find(total >1);
%        for i = 1:size(tri,1)
%            if hasOuterEdge
%                v1 = totalBoundary( tri(i,1));
%                v2 = totalBoundary( tri(i,2));
%                v3 = totalBoundary( tri(i,3));
%                
%                ismember( v1, side1,'rows')&& ismember( v1, side2,'rows')
%                          triangle(i).V1Side = 3;
%        
%        pt1 = 0;
%        for i = 1:size(outer,1)
%            if ismember(i,refPtsOnOuter)&& pt1 == 0
%                pt1 = i;
%            end
%            if ismember(i,refPtsonOuter)&& pt1 ~= 0
%                pt2 = i;
%            end
%        end
%        
%        %Look for a triangle in the mesh that contains point 1 and point 2
%        %and point 3 is on inner or opposite side of pt1 and pt2
%        for i = 1:size(refBoundaryPts
%            
%                
%                
%        
%        for i = 1:
       

%             %% Create triangle structure        
%              triangle = struct();
%              for i = 1:size(tri,1)
%                 v1 = totalBoundary(tri(i,1),:);
%                 v2 = totalBoundary(tri(i,2),:);
%                 v3 = totalBoundary(tri(i,3),:);
%                 triangle(i).Vertex1 = v1;
%                 triangle(i).Vertex2 = v2;
%                 triangle(i).Vertex3 = v3;
%                 triangle(i).Vref1 = tri(i,1);
%                 triangle(i).Vref2 = tri(i,2);
%                  triangle(i).Vref3 = tri(i,3);
%                 V12 = [[v2(1) - v1(1), v2(2) - v1(2)],0];
%                 V21 =  - V12;
%                 V13 = [[v3(1) - v1(1), v3(2) - v1(2)],0];
%                 V31= - V13;
%                 V23 = [[v3(1) - v2(1), v3(2) - v2(2)],0];
%                 V32 =  - V23;
%                 triangle(i).Vector12 = V12; 
%                 triangle(i).Vector21 = V21;
%                 triangle(i).Vector13 = V13;
%                 triangle(i).Vector31 = V31;
%                 triangle(i).Vector23 = V23;
%                 triangle(i).Vector32 = V32;
%                 triangle(i).Angle213 = (180/pi)*(atan2(norm(cross(V12,V13)),dot(V12,V13)));
%                 triangle(i).Angle123 =  (180/pi)*(atan2(norm(cross(V21,V23)),dot(V21,V23)));
%                 triangle(i).Angle132 =  (180/pi)*(atan2(norm(cross(V31,V32)),dot(V31,V32)));
%                 triangle(i).L12 = norm(triangle(i).Vector12);
%                 triangle(i).L13 = norm(triangle(i).Vector13);
%                 triangle(i).L23 = norm(triangle(i).Vector23);
% 
%                 %Determine the boundary that each vertex is on
% 
%                    %V1
%                    if ismember( v1, side1,'rows')&& ismember( v1, side2,'rows')
%                          triangle(i).V1Side = 3;
%                    elseif  ismember( v1, side1,'rows')
%                        triangle(i).V1Side = 1;
%                    elseif ismember( v1, side2,'rows')
%                        triangle(i).V1Side = 2;
%                    elseif ismember( v1, inner,'rows')
%                        triangle(i).V1Side = 0;
%                    else
%                           triangle(i).V1Side = 4;
%                    end
% 
%                    %V2
%                    if ismember( v2, side1,'rows')&& ismember( v2, side2,'rows')
%                          triangle(i).V2Side = 3;
%                    elseif  ismember( v2, side1,'rows')
%                        triangle(i).V2Side = 1;
%                     elseif ismember( v2, side2,'rows')
%                        triangle(i).V2Side = 2;
%                     elseif ismember( v2, inner,'rows')
%                        triangle(i).V2Side = 0;
%                     else
%                           triangle(i).V2Side = 4;
%                     end
% 
%                    %V3
%                    if ismember( v3, side1,'rows')&& ismember( v3, side2,'rows')
%                          triangle(i).V3Side = 3;
%                    elseif  ismember( v3, side1,'rows')
%                        triangle(i).V3Side = 1;
%                     elseif ismember( v3, side2,'rows')
%                        triangle(i).V3Side = 2;
%                     elseif ismember( v3, inner,'rows')
%                        triangle(i).V3Side = 0;
%                     else
%                         triangle(i).V3Side = 4;
%                     end
% 
%                     %Determine if the triangle is an external triangle
%                      if   (triangle(i).V1Side == triangle(i).V2Side...
%                            && triangle(i).V2Side == triangle(i).V3Side)...
%                        ||...
%                           (triangle(i).V1Side == triangle(i).V2Side...
%                            && triangle(i).V3Side == 3 )...
%                        ||...
%                            (triangle(i).V1Side == triangle(i).V3Side...
%                            && triangle(i).V2Side == 3 )...
%                         ||...
%                            (triangle(i).V2Side == triangle(i).V3Side...
%                            && triangle(i).V1Side == 3 )...
% 
%                           triangle(i).IsExternal = 1;
%                      else
%                          triangle(i).IsExternal = 0;
%                      end
% 
%                      %Determine if it is an end triangle
%                      if triangle(i).V1Side == 3 ...
%                          || triangle(i).V2Side == 3 ...
%                          || triangle(i).V3Side == 3
%                             triangle(i).IsEnd = 1;
%                      else
%                          triangle(i).IsEnd = 0;
%                      end
% 
%                      %Determine the edge along the boundary for right triangles
%                      if ~triangle(i).IsExternal && ~triangle(i).IsEnd
%                        if triangle(i).V1Side == triangle(i).V2Side
%                            triangle(i).IsRightTriangle = 1;
%                            triangle(i).BoundaryEdge = [1,2];
%                            if triangle(i).L13 < triangle(i).L23
%                               triangle(i).RightEdge = [1,3];
%                               a  = triangle(i).Vertex1;
%                               b  = triangle(i).Vertex3;
%                               triangle(i).Midpoint = [(a(1)+ b(1))/2,(a(2)+b(2))/2];
%                            else
%                               triangle(i).RightEdge = [2,3];
%                               a  = triangle(i).Vertex2;
%                               b  = triangle(i).Vertex3;
%                               triangle(i).Midpoint = [(a(1)+ b(1))/2,(a(2)+b(2))/2];
%                            end
%                        elseif triangle(i).V1Side == triangle(i).V3Side
%                            triangle(i).IsRightTriangle = 1;
%                            triangle(i).BoundaryEdge = [1,3];
%                            if triangle(i).L12 < triangle(i).L23
%                               triangle(i).RightEdge = [1,2];
%                               a  = triangle(i).Vertex1;
%                               b  = triangle(i).Vertex2;
%                               triangle(i).Midpoint = [(a(1)+ b(1))/2,(a(2)+b(2))/2];
%                            else
%                               triangle(i).RightEdge = [3,2];
%                               a  = triangle(i).Vertex3;
%                               b  = triangle(i).Vertex2;
%                               triangle(i).Midpoint = [(a(1)+ b(1))/2,(a(2)+b(2))/2];
%                             end
%                        elseif triangle(i).V2Side == triangle(i).V3Side
%                            triangle(i).BoundaryEdge = [2,3];
%                            triangle(i).IsRightTriangle = 1;
%                            if triangle(i).L12 < triangle(i).L13
%                               triangle(i).RightEdge = [2,1];
%                               a  = triangle(i).Vertex2;
%                               b  = triangle(i).Vertex1;
%                               triangle(i).Midpoint = [(a(1)+ b(1))/2,(a(2)+b(2))/2];
%                            else
%                               triangle(i).RightEdge = [3,1];
%                               a  = triangle(i).Vertex3;
%                               b  = triangle(i).Vertex1;
%                               triangle(i).Midpoint = [(a(1)+ b(1))/2,(a(2)+b(2))/2];
%                            end 
% 
%                        else
%                            triangle(i).BoundaryEdge = [0,0]; 
%                            triangle(i).IsRightTriangle = 0;
%                            triangle(i).Midpoint = [0,0];
%                        end
%                      else
%                        triangle(i).BoundaryEdge = [0,0]; 
%                        triangle(i).IsRightTriangle = 0;
%                        triangle(i).Midpoint = [0,0];
%                      end    
% 
% 
% 
%                 if triangle(i).IsExternal == 0 && triangle(i).IsEnd == 0 
% 
%                     if    triangle(i).V1Side == 4 ...
%                        || triangle(i).V2Side == 4 ...
%                        || triangle(i).V3Side == 4
%                             triangle(i).IsInternal = 1;
%                     elseif triangle(i).V1Side ~= triangle(i).V2Side ...
%                           && triangle(i).V2Side ~= triangle(i).V3Side ...
%                           && triangle(i).V1Side ~= triangle(i).V3Side
%                           triangle(i).IsInternal = 1;
%                     else 
%                         triangle(i).IsInternal = 0;
%                     end
%                 else
%                     triangle(i).IsInternal = 0;
%                 end
% 
%              end
% 
% 
% 
%                 disp('Boom');
% 
%             %Plot bw image
%              figure(1), imshow(bwImage) 
% 
%             %Plot inner an outer boundary together
%             figure(2)
%             plot(totalBoundary(:,1),totalBoundary(:,2));
%             
% 
%             %Plot the entire boundary and mesh
%             figure(3)
%             triplot(tri, totalBoundary(:,1), totalBoundary(:,2));
% 
%             %Plot inner and outer boundary with sides shown
%             figure(4)
%             plot(side1(:,1), side1(:,2), 'g');
%             hold on
%             plot(side2(:,1), side2(:,2), 'm');
%             hold on
%             plot(xInner, yInner, 'b'); 
%             hold off
%             
%            %Plot with external mesh removed  
%            remove = zeros(size(tri,1),1);
%            for i = 1:length(triangle)
%                if triangle(i).IsExternal 
%                    remove(i,1)  = 1;
%                end   
%            end
%            index = find(remove);
%            trimod = removerows(tri,'ind', index);
%            figure(5)
%            triplot(trimod, totalBoundary(:,1), totalBoundary(:,2));
% 
% 
%            %Plot with external meash and internal triangles removed
%             remove = zeros(size(tri,1),1);
%            for i = 1:length(triangle)
%                if triangle(i).IsExternal == 1 %|| triangle(i).IsInternal ==1
%                    remove(i,1)  = 1;
%                end   
% 
%            end
%            index = find(remove);
%            trimod = removerows(tri,'ind', index);
%            figure(6)
%            triplot(trimod, totalBoundary(:,1), totalBoundary(:,2));
% 
%             %Plot with external mesh, internal triangles, and triangles
%             %sharing vertices with internal triangles.
%             %Plot with external meash and internal triangles removed
%             remove = zeros(size(tri,1),1);
%             vertexRef = [];
%            for i = 1:length(triangle)
%                if triangle(i).IsExternal == 1
%                    remove(i,1)  = 1;
%                end
%                if triangle(i).IsInternal == 1  
%                    remove(i,1) = 1;
%                    vertexRef = [vertexRef,triangle(i).Vref1];
%                    vertexRef = [vertexRef,triangle(i).Vref2];
%                    vertexRef = [vertexRef,triangle(i).Vref3];
%                    for j = 1:size(tri,1)
%                         if ismember( tri(j,1), vertexRef) || ismember( tri(j,2), vertexRef) || ismember( tri(j,3), vertexRef) 
%                             remove(j,1) = 1;
%                         end
%                     end
% 
% 
%                end   
% 
%            end
% 
% %            %Estimate the skeleton in missing sections
% %             for i = 1:length(triangle)
% %                if triangle(i).IsExternal == 1
% %                    remove(i,1)  = 1;
% %                end
% %                if triangle(i).IsInternal == 1  
% %                   
% %                    vertexRef = [vertexRef,triangle(i).Vref1];
% %                    vertexRef = [vertexRef,triangle(i).Vref2];
% %                    vertexRef = [vertexRef,triangle(i).Vref3];
% %                    for j = 1:size(tri,1)
% %                         if ismember( tri(j,1), vertexRef) || ismember( tri(j,2), vertexRef) || ismember( tri(j,3), vertexRef) 
% %                          tNormal = triangle(i).RightEdge ;
% %                          tBoundary = triangle(i).BoundaryEdge;
% %                          [index, loc] = ismember(tBoundary,tNormal);
% %                          boundaryTerminal1 =  totalBoundary(:tNormal(loc));
% %                               
% %                                   tBoundary(1) a  = triangle(i).Vertex2;
% %                               b  = triangle(i).Vertex1;
% %                               triangle(i).Midpoint
% %                                
% %                         end
% %                     end
% % 
% % 
% %                end   
% 
%            end
%            index = find(remove);
%            trimod = removerows(tri,'ind', index);
%            figure(6)
%            triplot(trimod, totalBoundary(:,1), totalBoundary(:,2));
% 
% 
%            %Plot the midpoints of the right edges
%            seqnum = 0;
%            midpoints = [];
%            for i = 1:length(triangle)
%                if triangle(i).IsRightTriangle && ~isequal(triangle(i).Midpoint,[0,0])...
%                        && remove(i) ~=1
%                    seqnum = seqnum + 1;
%                    midpoints(seqnum,:) = triangle(i).Midpoint;
%                end
%            end
%            figure(7)
%            plot( midpoints(:,1), midpoints(:,2),'LineStyle','none','Marker','*')
% 
% 
% 
% %             %% Estimate of inner boundary
% %            if length(B) > 1;
% %                %Sample the outer boundary
% %                 for i = 1:5:size(xOuter,1)
% %                  xSampled = xOuter;
% %                  ySampled = yOuter;
% %                 end
% % 
% % 
% %                 make_plot = 0;
% %                 flag1 = 0;
% %                 d = 20;
% %                 [xInnerHat, yInnerHat, ~, ~, R, unv, concavity, overlap]=parallel_curve(xSampled, ySampled, d, make_plot, flag1);
% % 
% % 
% %                 triHat = delaunay([xOuter;xInnerHat], [yOuter;yInnerHat]); 
% % 
% %         %         x_inner= x_inner(1:4);
% %         %         y_inner = y_inner(1:4);
% %         %         xx= [x_inner(1):-1:x_inner(4)];
% %         %         yy = spline(x_inner,y_inner,[x_inner(1):-1:x_inner(4)]);
% % 
% %                 figure(8)
% %                 plot( midpoints(:,1), midpoints(:,2),'LineStyle','none','Marker','*')
% %                 hold on;
% %                 triplot(triHat, [xOuter;xInnerHat], [yOuter;yInnerHat]);
% %                 hold off
% %            end
% 
%              %legend({'Curve', 'Inner Parallel', }, 'location', 'Best');
% 
%             
%         
%         
%        
        pause;
    end


